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The Turing mechanism 1 for the production of a broken spatial symmetry in an initially 
homogeneous system of reacting and diffusing substances has attracted much interest as a 
potential model for certain aspects of morphogenesis 2 ^ such as pre-patterning in the embryo, 
and has also served as a model for self- organization in more generic systems. 5 The two 
features necessary for the formation of Turing patterns are short-range autocatalysis and 
long-range inhibition 6 ' 7 which usually only occur when the diffusion rate of the inhibitor is 
significantly greater than that of the activator. This observation has sometimes been used 



to cast doubt on applicability of the Turing mechanism to cellular patterning since many 
messenger molecules that diffuse between cells do so at more-or-less similar rates. Here we 
show that stationary, symmetry-breaking Turing patterns can form in physiologically realistic 
systems even when the extracellular diffusion coefficients are equal; the kinetic properties of 
the "receiver" and "transmitter" proteins responsible for signal transduction will be primary 
factors governing this process. 

The class of mechanisms we study is schematized in Fig. 1. Below we give a specific 
example for the purpose of illustrating the generic features of the whole class of such systems. 
The model is realistic in the sense that it makes use of reactions and reaction kinetics of 
a type which can be found in every living cell. In the model, signal transduction across 
the cell membrane connects the genetically controlled biochemical reactions in the cytosol 
to the production of messenger molecules in the extracellular matrix. These then diffuse at 
approximately equal rates, thereby coupling the reactions taking place inside cells at different 
points in space. We have found that this coupling can allow Turing patterns to spontaneously 
form in collections of cells even if the messenger diffusion rates (D) are identical (or very 
similar). We therefore believe the chemically realistic class of mechanisms we describe here 
could be of some relevance to understanding certain aspects of cellular morphogenesis. 

The simplified set of reactions taking place in the cytosol involve an "activator" substance 
A and an "inhibitor" substance / that are synthesized by the cell at the rates ca and C/ 
respectively, and in turn are broken down by the cell at the rates Xa and A/ as shown in Fig. 



1, where 

d[A]/dt = c A - X A [A], d[I]/dt = a - A,[/]. (1) 

Activator is also produced by the autocatalytic process 

A + E x ^ E X A (K x = [E X ][A\/[E X A\), t/> + E X A _^ E X A + A, (2) 

d[A]/dt = V X [A]/(K X + [A]), V x = [E x ] c k x , 

where we will use ip as a shorthand for substances from the large and constant concentration 
of assorted biochemical building blocks present in the cell that are used to synthesize A and 
/. The activator also catalyzes the production of more inhibitor through the set of reactions 

k 2 
A + E 2 ^ E 2 A (K 2 = [E 2 }[A]/[E 2 A\), t/> + E 2 A _^ E 2 A + /, (3) 

d[I]/dt~V 2 [A]/K 2 , V 2 = [E 2 ] c k 2 , 

where for simplicity we have assumed that this reactions operates well below saturation, 
with [A] <C K 2 . Lastly, the inhibitor suppresses the activator through the set of reactions 

A + E 3 # E 3 A (K 3 = [E 3 ][A]/[E 3 A]), (4) 

h 

I + E 3 A _^ E 3 AI -+ E 3 ^I -> E 3 + i) + /, 



-d[A]/dt = V 3 [A][I]/(K 3 + [A]), V 3 = [E z ] c k z - 

The rates of the catalytic process are given by the Michaelis-Menten 8 kinetics, where Ki 
are the Michaelis constants, ki the catalytic rates and V* the limiting rates for the catalytic 
processes. The E^s are catalysts of the reactions where [Ei] are the total concentrations 
of these catalysts produced by the genetic machinery of the cell. Because of the need for 
enzymes, we assume these reactions, including the breakdown of A and I, can only occur 
inside the cell. 

In addition to being biologically reasonable, our model differs from the classical Turing 
model in two important respects. First, it takes into account the fact that genetically regu- 
lated biochemical processes will take place inside the cell, while cell-to-cell interaction must 
involve some form of trans-membrane signalling, since most biomolecules are highly insolu- 
ble in the lipid matrix of the cell membrane. Secondly, the diffusion of the two messenger 
molecules through the extracellular matrix take place at identical rates D. 

Since pre-patterning is likely to involve switching on different sets of genes in different 
cells, some mechanism must exist for signalling between the genes. Our model makes use 
of signal transduction kinetics diagrammed in Fig. 1, where a a and olj are the rates of 
production of the messenger molecules Ma and Mj. This process is mediated by transmit- 
ter proteins Ta and Tj in the membrane. In the figure, signal transduction is viewed as 
the transformation of A and I into Ma and M/ respectively (though transformation into 
messenger molecules is not a requirement of the model). Conversely, the Pa and j3j are the 
rates of transformation of Ma and M r into A and /, a process mediated by Ra and Ri, the 



receiver proteins. This simplified picture can stand for more complex signal transduction 
mechanisms since the basic features responsible for the self-organization of the system are 
generic. For example, the signal transduction could involve the binding of A or 7 to specific 
receptor proteins on the intracellular side of the membrane, triggering the release of messen- 
ger molecules on the extracellular side of the membrane which then carry the signal to other 
cells. 

The cell essentially has nearly complete freedom to control the signal transduction ki- 
netics, whose associated rates could quite easily vary by many orders of magnitude. This is 
to be strongly contrasted with the situation where the diffusion coefficients themselves are 
required to differ by a large amount - the standard route to Turing patterns - since it is 
difficult to imagine realistic situations where this holds true. 

The biological mechanism we propose here is analogous to the mechanism proposed as 
an explanation for the chlorite-iodide-malonic acid-starch reaction system, which was the 
first chemical system in which Turing patterns were observed experimentally. 9 Lengyel and 
Epstein 10 ' 11 showed that the existence of Turing patterns in this system, despite the similar 
diffusion coefficients of the reactants, is a consequence of the binding and unbinding of 
iodine to starch molecules that have been immobilized by the gel which serves as an "inert" 
medium of the reaction. What we are proposing here is a biologically realistic and necessary 
mechanism that is able to provide just this type of effect in living cells. 

Setting A(x,t) = [A], 7(x,£) = [7] M A (x,t) = [M ] and M 7 (x,t) = [M,] to be the 
concentration at a given time t and position x, we can write down the following set of four 



reaction-diffusion equations 



ViA V 3 AI .,.,,, 



d t A = c A - X A A + — ±— - - — i— - - a A A + /5 A M A - eVM, (5) 

9 t 7 = C/ - A/7 + %^ - a// + /3/M/ - eV 2 /, (6) 

9 t M A = a^ - /^M 7 - DV 2 M A , (7) 

d t M 7 = ajl - fijMi - DV 2 Mj. (8) 

Here e is a very small (or negligible) rate of diffusion arising from "leakage" though the cell 
membrane which occurs very rarely. We introduce this small quantity for a mathematical 
reasons in order to give the first two equations a spatial scale and to remove the singular 
nature of the limit where e = 0. However we emphasize that the model still forms patterns 
even when e = (see Fig. 2). 

For the purposes of keeping our illustration as simple as possible we will specifically 
consider only the case where c A = cj = 1 , a A = Pa = 1 and olj = fii = p. We further 
set limiting velocities of the catalyzed reactions to V\ = 500, V 2 = 1, and V 3 = 60, and the 
parameters A^ = 1/100, A/ = 7, K\ = 100, K 2 = 1, and K 3 = 1/10. For this model the 
stationary homogenous state is A(x.) = M A (x.) ~ 2.46 and 7(x) = M/(x) ~ 2.25. A standard 
linear stability 2 analysis of Eqs. (1-4) about this state tells us that an initial perturbation of 
the uniform state with a given wave vector k will have an amplitude that grows (or shrinks) 



in proportion to the factor e * , where A& are eigenvalues of the stability matrix 
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and can be real or come in complex conjugate pairs. Each value of k corresponds to a spatial 
mode with wavelength L = 2ny / D/k. The real parts of the eigenvalues corresponding to 
this mode, 9?(Afc), determine the stability of perturbations away from the homogeneous state 
with the corresponding wavelength: 5R(A^) < indicates the mode is stable, and 5R(A^) > 
indicates an unstable mode which will grow when a small random fluctuation displaces it 
from equilibrium. Thus the existence of a wave number k such that A& > means patterns 
will form. Turing patterns occur when the largest eigenvalue is real (S(A&) = 0) which is 
the case for the parameters we have chosen here. The largest eigenvalue has a negative real 
part (9?(Afc) < 0) when D — > showing that the system is globally stable in the absence of 
diffusion. 

The type of patterns that form will depend to a large degree on the spectrum of unstable 
k modes; it especially depends on the mode with the greatest A&, but also on the range of 
other unstable modes. As was mentioned above, we included a small diffusion effect through 
the membrane with diffusion coefficient e for mathematical reasons in order to remove the 
singular nature of the equations when the diffusion vanishes. In Fig. 2(B) we show the 
effect of varying e on the frequency spectrum. As e — >■ the spectrum of unstable modes is 
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broadened towards the higher frequency modes until at e = 0, the maximum k value in the 
spectrum goes to infinity. Physically this means that shorter and shorter wavelengths will 
come to dominate, and at some point the wavelength of the patterning will be limited by 
the finite dimension of the cells. 

Fig. 3 shows some of the ordered patterns that emerge spontaneously in the system 
we studied for varying values of p ranging from (A) "honeycomb" (p = 6) to (B) stripes 
(p = 11.5) to (C) spots (p = 100). Our system is capable of supporting more or less the same 
patterns as an ordinary Turing mechanism that uses widely differing diffusion coefficients. 
The principal difference is that our mechanism provides a physiologically and chemically 
feasible route by which a wide range of patterns could arise. 

Previous work has shown that pattern formation can occur in some systems with equal 
diffusion coefficients, but only if there is some kind of initial asymmetry such as a finite 
size perturbation 12 or an external advective flow. 13 ' 14 These mechanisms might provide 
useful models for subsequent morphogenetic events in the embryo where some spatial genetic 
patterning or chemical gradients have already been set up, but still does not provide an 
adequate explanation of initial pre-patterning events which may start from a completely 
homogeneous initial state. The existence of an organizing center begs the question of how 
the center itself formed. 

Our model provides a possible, physiologically realistic route to symmetry breaking in- 
stabilities in cellular systems, and we hope that it will provide a useful context in which 
to explore the possible relevance of Turing mechanisms to cellular morphogenesis. For in- 



stance, our model makes a clear distinction between reacting and diffusing entities and thus 
provided somewhat different expectations for making experimental observations. The mes- 
sengers could be almost any type of molecule, and are not required to have any complex 
reactive chemical properties. Rather, it is the properties of the transmitter and receiver 
proteins in the signal transduction pathway that will control the pattern formation. As a 
consequence it may be unproductive to search for specific "morphogens" since the biochem- 
ical substances that differentiate spatial patterns in a collection of cells may be unrelated to 
the substances that actually mediate the cell-to-cell communication. 

M.M. is supported by a Whitaker Foundation award. E.R. is supported by a National 
Science Foundation Graduate Fellowship. 
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Figure Legends 

Figure 1: General kinetic schematic for model. The three basic kinetic elements of 
our model are shown. The simplified reactions in the cytosol involve an activator substance 
A and inactivator i" which are created at constant rates by the cell and likewise broken down 
at the rates A^ and A/. In this case trans-membrane signal transduction takes the form of 
transformation of A and / in the cytosol (mediated by the membrane proteins T A and Tj) 
into corresponding messenger molecules Ma and Mi in the extracellular matrix. Likewise 
Ra and R r mediate the reverse transformation. Both of the messenger substance diffuse at 
rate D through the extracellular matrix. 

Figure 2: Linear stability spectrum of the model. The dotted lines indicate regions 
where the eigenvalue with the largest real part is complex. Parametet given in the text. (A) 
shows the spectrum of stability eigenvalues for several values of the critical parameter p. For 
p < 3.5 all the eigenvalues have negative real parts, and the homogeneous distribution of 
reactants is stable. Turing patterns set in when p > 3.5, and a range of frequency modes 
become unstable, giving rise to spontaneous ordering. (B) shows the effect on the spectrum 
when p = 5 as e is varied. As e — >■ 0, the most unstable mode shifts to higher and higher 
frequencies, but ordering will still take place. 

Figure 3: Turing Patterns (A) Concentrations of activator and inhibitor within cells at 
t = 125 for p = 6. Concentration were set to zero at t = and allowed to build up naturally, 
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the situation that would arise if the genes for the production of A and / were suddenly 
switched on at t = 0. The activator is shown as shades of green (black indicating and green 
the maximum value, here 32.0), and the inhibitor as shades of red (black = to red = 21.3) as 
shown in the legend (D). The two are superimposed so that yellow indicates the presence of 
both activator and inhibitor. D = 600, e = D/100; the grid size is 100. (B) Concentrations of 
A and / at t = 40 for p = 11.5 and other parameters as above. A max = 53.8 and I max = 28.4. 
(C) Concentrations of A and / at t = 40 for p = 100. A max = 117.6 and I max = 35.3. 
Numerical Method: The reaction component of the equations was integrated using a 
forward-Euler method. The Laplacian (diffusion) term was implemented using the following 
conservative method: for each pair of locations {i, j}, Ai(t + 1) = Ai(t) + AA; Aj(t + 1) = 
Aj(t)-AA where AA = M ij D(A i (t+l)f(i,t)-A j (t+l)f(j,t))At and M tj is the connectivity 
matrix. The noise function / is a uniformly distributed random variable, ranging from 
to 10~ 5 and models the effect of very weak fluctuations of the type that are required to 
initiate the initial instability. We used a von Neumann neighborhood and toroidal boundary 
conditions. 
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